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Abstract In this paper we study a history matching approach that consists of find- 
ing stable approximations to the problem of minimizing the weighted least-squares 
functional that penalizes the misfit between the reservoir model predictions G(u) 
and noisy observations y. In other words, we are interested in computing ifl = 
argmin i(G x _1 / 2 (y — where F is the measurements error covariance, Y 

is the observation space and X is a set of admissible parameters. This is an ill-posed 
nonlinear inverse problem that we address by means of the regularizing Levenberg- 
Marquardt scheme developed in [7 8]. Under certain conditions on G, the theory 
of C][8 1 ensures convergence of the scheme to stable approximations to the inverse 
problem. We propose an implementation of the regularizing Levenberg-Marquardt 
scheme that enforces prior knowledge of the geologic properties. In particular, the 
prior mean u is incorporated in the initial guess of the algorithm and the prior er- 
ror covariance C is enforced through the definition of the parameter space X. Our 
main goal is to numerically show that the proposed implementation of the regulariz- 
ing Levenberg-Marquardt scheme of Hanke is a robust method capable of providing 
accurate estimates of the geologic properties for small noise measurements. In ad- 
dition, we provide numerical evidence of the convergence and regularizing results 
predicted by the theory of 00 for a prototypical oil-water reservoir model. The 
performance for recovering the true permeability with the regularizing Levenberg- 
Marquardt scheme is compared against the more standard techniques for history 
matching proposed in [16,21,22,19]. Our numerical experiments suggest that the 
history matching approach based on iterative regularization is robust and could po- 
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tentially be used to improve further on various methodologies already proposed as 
effective tools for history matching in petroleum reservoirs 

1 Introduction 

History matching is the process of modifying parameters (inputs) of a reservoir model 
so that the model response (output) matches production data. The adjusted parameters 
during the history matching process are often geologic properties of the subsurface 
whose lack of information gives rise to uncertainty in the predictions of the reservoir 
model. The parameters obtained by means of history matching are aimed to provide 
better predictions of the reservoir performance; these can potentially be used for opti- 
mal reservoir management and monitoring of the reservoir. Given the potential impact 
of the history matching process in the optimal production of hydrocarbons, numer- 
ous techniques have been proposed and widely investigated in the last decades. For 
a recent review of history matching techniques we refer the reader to |20| . Standard 
approaches are presented in detail in the monograph 1 19 1. 

Let us denote by u the unknown parameter (geologic properties) in a reservoir 
whose dynamics are described with a parameter-to-output operator G : X — > Y that 
maps the space of admissible parameters X to the observation space Y. In this paper 
we study the history matching problem posed by the minimization of the weighted 
data misfit defined by 

®( u )= l -\\r-^(y*-G{u))\\l (1) 

where is the production data and F is the measurements error covariance. For 
reservoir modeling applications, the evaluation of the forward operator G(u) involves 
the solution of highly nonlinear system of PDEs whose differential operators have 
spatially varying coefficients related to the geologic parameter u. Therefore, the op- 
erator G is typically compact which from standard theory |5| implies that the mini- 
mization of ([T]i is ill-posed in the sense of stability. In other words, a small perturba- 
tion of may correspond to large deviations from the corresponding solutions to the 
minimizer of ([l} ll5l [T5l . The aforementioned lack of stability can lead to the diver- 
gence of standard optimization schemes used to solve the least-squares problem ([T| 
ifTTl Hl. In this paper we propose a computational approach for the solution of ([TJ by 
means of the regularizing Levenberg-Marquardt (LM) scheme developed by Hanke 
in El|8 |. The regularizing LM scheme belongs to the class of so-called iterative reg- 
ularization techniques designed to compute stable approximation to inverse ill-posed 
problems like the one posed by the minimization of In contrast to the other ap- 
proaches where the problem is first regularized (e.g. by Tikhonov's method) and then 
optimized with a standard solver, with iterative regularization techniques, the aim is 
to regularize the problem within an algorithm that also provides an approximation 
to a minimizer of ([TJ. In other words, iterative regularization schemes provide stable 
estimates that converge to a minimizer of <P in the limit of small noise. A review of 
the analysis and applications of iterative regularization techniques can be found in 
fl4l . For the regularizing LM scheme, the mathematical analysis of the convergence 
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and regularizing properties is developed in [7 8|. This mathematical framework of 
the LM scheme motivates the implementation that we propose for history matching 
in petroleum reservoirs. The main objective of this paper is to numerically show that 
the proposed implementation of the regularizing LM scheme is a robust methodology 
for generating accurate estimates of geologic parameters given production data with 
small noise. Furthermore, we provide numerical comparisons of the performance of 
the proposed implementation with respect to one of the most standard approaches 
for deterministic history matching. While the regularizing LM scheme aims at solv- 
ing a history matching problem posed differently from standard methods, there exist 
some technical similarities between our implementation and the standard approaches. 
We exploit those similarities to provide guidelines for a straightforward implementa- 
tion of the proposed algorithm given routines and codes from standard optimization 
methods 

The paper is organized as follows. Relevant literature is discussed in Section[2] In 
Section[3]we introduce the application of the regularizing LM scheme of [78| to the 
history matching problem. The fundamental theoretical aspects of the regularizing 



LM scheme are discussed in subsection 3.1 Computational aspects relevant to the 



implementation of the regularizing LM scheme are presented in subsection 3.2 In 
subsection 3.3 we discuss one of the standard approaches for history matching based 
on the method used in [ 16,21 ,22, 19]. In subsection |3.4| we show fundamental simi- 
larities between the proposed implementation of the regularizing LM scheme and the 
aforementioned standard approach. In Section[4]we display numerical experiments to 
show the capabilities of the regularizing LM scheme for generating stable approxima- 
tions to the proposed history matching approach. Our implementation of the regular- 
izing LM scheme is applied to an incompressible oil-water reservoir model described 
in the Appendix. The regularizing properties of the LM scheme with respect to the 
noise level are studied in subsection 4.2 The performance of the LM scheme with 
respect to relevant tunable parameters is illustrated in subsection 4.3 The efficacy 
of the regularizing LM scheme for different choices of variance in the prior covari- 
ance operator is investigated in Section 4.4 Finally, the proposed regularizing LM 
scheme is compared against the standard method of 1 16,21 22 19 1. Conclusions and 
final remarks are provided in Section[5] 



2 Literature Review 

For history matching applications posed in terms of ([TJ, regularization has been typ- 
ically addressed by reparameterizing the geologic properties with a small number of 
parameters (see [20], section 3.2 and references therein). A parameterization in terms 
of a finite dimensional (hence compact) set, ensures the well-posedness of ([TJ [13|. 
However, for history matching applications, only problems parameterized with very 
few (of the order of 10) parameters have been treated by minimizing a functional 
like ([TJ. For reservoirs with highly heterogeneous geologic properties, thousands or 
even millions of parameters are required to fully resolve relevant geologic features. 
For those reservoirs, parameterizing the geologic properties with a small number of 
parameters may not be possible. In that case, minimizing ([TJ with standard optimiza- 
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tion techniques may diverge due to the lack of stability described above. One of the 
most standard approaches for history matching that addresses the ill-posedness of the 
inverse problem is to minimize |20| 

J{u)^\\\r- l l 2 (y^-G{u))\\l + \\\C- l l\u-ml (2) 

which can be thought as Tikhonov regularization of Q [ 19]. Indeed, the second term 
in the right hand side of (|2| 

R{u)^ l -\\C- l l 2 {u-u)\\ 2 x . (3) 

is a regularization term that alleviates the ill-posedness of the inverse problem ([lj. 
Under some assumptions on G, the theory of Tikhonov regularization for nonlinear 
problems ensures that (|2| has a solution which is continuous with respect to J5] 
Theorem 10.2]. Therefore, the problem is well-posed and any standard optimization 
technique can be implemented for solving (j2j). However, as the size of R(u) decreases, 
the regularization properties of (|2jl rely on the proper size of R{u) relative to 'P(u) 
J5] Theorem 10.3]. Intuitively, if the regularization term R(u) is "too small" (under- 
estimated) compare to <P(u), the regularization provided by R(u) may not suffice, 
resulting in the lack of stability. On the other hand, if R(u) is "too large" (overesti- 
mated), the minimizer Q may produce estimates that lackfidelity due to a potential 
poor data match. For a given fixed error covariance matrix F, then the relative size of 
R(u) with respect to <P(u) is determined by the covariance operator C. For reservoir 
applications the standard practice is to select C based on geologic data which has no 
connection to the stability and fidelity issues of the inverse problem described above. 
Therefore, for some choice of C, the potential risk of instabilities and lack of fidelity 
in minimizing |2} may arise. History matching applications where the minimization 
of |2) with the Gauss-Newton led to instabilities have been reported in B16II19I . In or- 
der to alleviate these instabilities, a standard Levenberg-Marquardt method has been 
proposed in [ 1 6 2 1 22 1 9 ] . The well-posedness of the minimization of Q is a funda- 
mental assumption for the application of those standard techniques. However, as we 
mentioned above, for some choices based on geological information of C may result 
in an insufficient regularization of the term R(u) which in turn, may lead to numerical 
instabilities in more general settings. 

By minimizing ([T]l instead of the implementation of the regularizing LM 
scheme that we propose in this paper avoids the potential lack of fidelity of the stan- 
dard approach previously discussed. While in the standard approach of minimizing 
(|2| the term R{u) enforces the prior geological knowledge, in our implementation of 
the regularizing LM scheme, the prior mean is incorporated in the initial guess of 
the iterative algorithm and the geological constraint imposed by the prior covariance 
is enforced in the definition of the parameter space X. It is fundamental to empha- 
size that the regularizing LM scheme applied to the minimization of ([T]) is a regu- 
larization technique that aims at producing stable computational approximations to 
a minimizer of ([TJ. Therefore, in contrast to the standard approach where the mini- 
mization of |2| is assumed well-posed and so standard optimization techniques can 
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be applied, the regularizing LM approach of Hanke postulates an algorithm that alle- 
viates the ill-posedness in the minimization of ([TJ while computing an approximation 
that converges to a minimizer of ([T]i for small observational noise (71IH1- 

Iterative regularization techniques such as the regularizing LM scheme, have been 
successfully used for the solution of a wide class of inverse problems in several dis- 
ciplines. In particular, for the inversion of data in subsurface flow models, in the 
authors study a simplified version of the regularizing LM scheme (see discussion of 
Section[3]) for the inversion of pressure in single-phase Darcy flow. In lfl2l a truncated 
regularizing Newton-Conjugate gradient is implemented to invert combined surface 
deformation and pressure data in a coupled flow-geomechanics problem. In this con- 
text of data inversion, the present work aims at extending the treatment in J9[ by 
using a two-phase (oil-water) reservoir model which, in contrast to the model stud- 
ied in [9|, is nonlinear with respect to the state variables (pressures and saturations). 
The increase in nonlinearity of the forward operator imposes severe challenges on the 
regularization properties of the technique under consideration. However, the present 
work offers the numerical evidence that the theory of Hanke may be applied to the 
forward operator that arises from the prototypical oil-water reservoir model that we 
consider in our numerical experiments. Further investigations of the regularizing LM 
scheme and other iterative regularization technique may lead to the development of 
efficient tools for history matching applications. 



3 Iterative regularization for history matching 

In this section we present the application of the regularizing LM scheme of HO 
for history matching by means of minimizing <£ defined in ([T| where G is an arbi- 
trary reservoir model that captures the flow dynamics perfectly. For our experiments 
of Section |4j we use a forward operator G that we derived from a prototypical in- 
compressible oil-water model (see Appendix). As we indicated in Section [T] due to 
the ill-posedness of the minimization of (JlJ, a regularization algorithm is required 
to compute stable solutions to the inverse problem. While a broad spectrum of iter- 
ative regularization techniques can be used, here we consider the regularizing LM 
technique because of the computational similarities with standard LM methods that 
are typically used for standard approaches in history matching that we described in 
Section [33] 

Assume we are provided measurements possibly corrupted by noise y 7 ! with the 
noise level denoted by tj and defined by 

where u' denotes the true geologic properties of the reservoir. Note that if we knew 
the truth u T , in the absence of observational error (i.e. tj — > 0), we would measure the 
model predictions of the truth y = G(u^). Therefore, tj in (j^j) is an upper bound for the 
observational noise. In practice, tj can be defined from the measurement information 
also used for defining the measurement error covariance F. 
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3.1 The regularizing Levenberg-Marquardt method 

The aim of the regularizing LM scheme is to compute stable approximations to a 
minimizer of ([TJ). In other words, we want to compute u-q such that — > u as rj — >• 0, 
where u is a minimum of ([TJ in the limit tj — > 0. The approximation u 1 ! is the limit 
of a finite sequence of estimates {«m}m=i computed as we now describe. Given, the 
estimate u} n at the mth iteration of the scheme, the aim is to construct an update 
u m+i = M '« +^ M >«' where the increment Au} n is obtained by solving 

y*-G{ul)=DG{ul)Aul, (5) 

where DG(u m ) is the Frechet derivative of G at u m . Note that ([5]> is a linearized version 
of the equation — G(u) satisfied by a minimizer u in the case that the minimum in 
|l]l corresponds to 0(u) = 0. Note that the truth increment defined by Au^ — — ii} n 
satisfies 

G{u f )-G{ul)-R(ul,iS)=DG(ul)Ai4 n , (6) 

where R(um,u') is the Taylor remainder of G at u r around u} n . As we mentioned 
in Section [TJ the ill-posedness in the minimization of (JJjl can be attributed to the 
compactness of the forward operator which is often encountered in PDE-constrained 
inverse problems lfT3~l . From the compactness of the Frechet derivative of a compact 
operator (0, it follows that the linear operator DG{u}n) is compact for each m 6 N. 
Therefore, the linear inverse problem |5]l also requires regularization. In the regulariz- 
ing LM scheme of Hanke [7 8], Tikonov regularization is applied to (|5]l by computing 

Aul(a) = argrriin wGX /f M (w>, a) (7) 

where 



(8) 



The choice of a is fundamental to ensure the proper regularization of the inverse 
problem. Hanke proposes a such that 

||r- 1 / 2 (/-G(4)-DG(4)4^(a))||2>p 2 ||r-V2(^_G(4))||2 (9) 

for some p € (0, 1) (note that Au^,(a) in ^ depends on a via (|tJ-(|8J). 

To gain further insight of the regularizing LM scheme as well as the selection of 
a, let us define 

d"» = y" - G{ul), d m = G{u^) - G{ul)-R{ul,u% 

DG{ul)=g m , u m = u m -u. (10) 

which applied to |5]l and |6| yields 

d^ m = g m Aul, d m =g m Aul (11) 



From definitions ( 10 1, expressions ([8]l and (|9]) become 



J LM 



;w,a)^i||r- 1 /2(^-_^4(a))||2 + Ia||C- 1 / 2 44(a)||2 (12) 
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and 

\\r-y\d^ m -g m Aul{a))\\\>p z \\r- l l 2 d^ m \\\ (13) 

respectively. Therefore, each iteration of the proposed scheme can be viewed as a 
Tikhonov regularization for the linear inverse problem of find Au} n given data d r,,m , 
where the latter is a noisy version of d m . Note that, from ( flO| ) it follows that 

\\r-'l 2 (yn -y-R{ul^))\\ Y = \\r- l l 2 (d m -d^)\\y (14) 

The regularizing LM scheme assumes that it is possible to find p € (0, 1) such that 

l|r- 1 / 2 (/- 3 '-^(«l« t ))||r = ||r- 1 / 2 K- £ /''' m )||F<p||r- 1 / 2 £ /''' m || F .(i5) 

The inequality in the previous expression implies that the size of the error in the data 
d m must be smaller than the size of the observations d^' m . It is certainly hopeless to 



invert data whose error is of the order of the size of the observations. The p in ( 15 i is 
used in expression ([9]) for choosing the regularization parameter a. Moreover, from 
(JT3J it is easy to see that the selection of a according to ( pj) implies 

\\r- l / 2 (d^ m - g m Aul(a))\\ 2 Y > \\r- l i 2 {d m -d^ m )\\y (16) 

which is the discrepancy principle applied to the inverse problem d^ m = g m Au} n . The 
discrepan cy p rinciple states that the estimate Au^ n {a) of the solution to the inverse 
problem ( 11 1 cannot produce an output g m Aii} n {a) whose associated error is better 
than the noise level. For a discussion of the discrepancy principle in the context of 
linear inverse problems the reader is referred to (6). Let us now denote by a m a 
solution to inequality ( fT~3] >. Then, the update of the regularizing LM scheme is defined 

by 



A +Aul[a m ) =ul + aiffmn weX J^f(w,a m ) (17) 



which provides a new estimate of the geologic properties. The existence of <x m is 
proven in |6, 14 1 (see also discussion below). The minimizer of (8 1 with a m given by 
|9]l provides a regularized solution to the linear inverse problem (5 I. Furthermore, the 
regularizing LM scheme is terminated provided the (k+ l)th iteration produces an 
estimate j such that 

\\r^ 1/2 (y 71 ~G(ul +l ))\\ Y < ttj < \\r~ 1/2 (y n -G{i$))\\y (18) 
for T > 1/p. The resulting estimate if 1 = u) + , is the desired stable approximation to 



the inverse problem of minimizing ([TJ. Expression ( 18 i is also an application of the 



discrepancy principle which, in this context, states that the data misfit obtained with 
approximation to the inverse problem ifl should not be smaller than the noise level 
77. Intuitively, if p sw 1 (with p < 1) then z can be chosen T ~ 1 which, in turn, may 
result in estimates that provide a good data fit. The regularizing LM scheme is now 
summarized below: 

Algorithm 1 (Regularizing Levenberg-Marquardt Scheme) Consider the initial es- 
timate Mq = u. Choose parameters p < 1 and T > 1/p. For each m= 1 , . . . , k, 
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(1) Forward simulation. Given uU t simulate the model response G(um). 

(2) Stopping rule (Discrepancy Principle). If (|7iS| holds then stop (i.e. m = k+1). 
Output: u} n . 

(3) Update. Define u^ +l according to (17) with J™ M defined in (jsjl and a" 1 is chosen 



according to the ((PJ. 

Remark 1 In the regularizing LM scheme, prior knowledge u of the unknown is in- 
corporated as the initial guess of the LM algorithm. In addition, the prior covariance 
C is included in the definition of the parameter space, which formally, can be defined 
as the completion of the original space X under the norm | \C~ l l 2 -\\x- This choice of 
the space is reflected in the second term of the right hand side of |8]l 

The application of the discrepancy principle for the selection of a in d9) as well as 



the termination of the algorithm ( 18 i are key aspects for the regularization properties 
of the regularizing LM scheme. In particular, we recall the following result proven in 
Theorem 2.3]. 

Theorem 1 (Hanke [7|) Let p £ (0,1) and T > 1/p. Assume that DG is locally 
bounded and that G satisfies 

\\G{u)-G(u)-DG(u)(u-u)\\ Y <C\\u~u\\x\\G(u)-G{u)\\ Y (19) 

locally in X. If uq is sufficiently close to a solution u* of G(w) = G(u*), then, the 
discrepancy principle (18 \ terminates the LM algorithm with parameters a from ([£]) 



after a finite number of iterations k(r\). Moreover, the corresponding approximations 
u k(ri) converge to a solution of G(u*) = G(u) as tj — > 0. 

Remark 2 From Theorem[T]we see that as 77 — >• 0, then the solution u* 1 computed with 
the regularizing LM scheme converges to u that satisfies G(m^) = G(u). Therefore, 
since from (j4j) y 1 ! — > G(u r ) as 77 — > 0, it then follows that u satisfies <P(u) = and so 
u 7 ! converges to a minimizer of <P in the limit of 77 — » 0. 

In (9) we implemented a particular case of Algorithm [TJ for the estimation of ab- 
solute permeability with a single-phase (linear) reservoir model. Instead of choosing 
a" as in (|9j, in the scheme of J9) the Tikhonov parameter was chosen constant a = 1 . 
This selection was sufficient to prove convergence in the same sense of Theorem [Tj 
For the work reported in this paper, we initially implemented the algorithm of [9] 
for the estimation of absolute permeability with the reservoir model of the Appendix. 
However, the need for an adaptive selection of a arose due to the highly nonlinear 
structure of the present forward model. While the rigorous application of Theorem[TJ 
for the forward operator G of the Appendix remains an open problem, our numerical 
results give evidence that confirms the regularizing properties predicted by Hanke's 
theory. 



3.2 Computational Implementation of the regularizing LM scheme 



In this section we discuss computational aspects of the regularizing LM scheme. Our 
main goal is to provide a reproducible computationally efficient algorithm for history 
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matching. We first notice that, for a fixed, the Euler-Lagrange equation associated to 
the minimization of ( 12 1 yields 

Aul{a) = \DG\ul)r- l DG{ul) + aC- l Y X DG*{ul)r- l {y^ -G{ul)). (20) 



where DG*(um) is the adjoint operator of DG(um). Expression (20i involves the in- 
version of the operator DG* (u m )r~ l DG(u m ) + aC~' in the space X. However, for 
the reservoir application under consideration, the dimension of the parameter space 
X is typically much larger than the dimension of the observation space Y. Therefore, 
for computational efficiency we consider the equivalence between ( 20 1 and 

Aul{a) =CDG*(ul)\DG(ul)CDG*(ul) + ar] (y^-G{ul)). (21) 

which in finite dimensions can be shown from the matrix lemmas of |fl9l Section 
7.4]. In the infinite-dimensional case, the equivalence between pO])-pT[) is only for- 
mal. Note that, assuming that the sensitivities DG(um) and DG*(i(m) are available, 



then either (20 1 or pi) can be easily computed for any given a. It is therefore clear 
that the computation of a in |9]) represents the main new aspect of the proposed im- 
plementation. However, the computation of a is fairly simple as we describe below. 
Let us define 

K^a)^\\r-^(y^-G(u^-DG(ul)AM<x))\\x (22) 
We substitute expression pT) in ([22]) and from simple computations it follows that 

K^a)^a 2 \\r^[DG(u^CDG^ u ^ + ar}-^-G(um 2 r (23) 

From this expression we find that fCm(a) is a continuous increasing function of a. 
Moreover, it can be shown |[T4l Chapter 4] that 



.\\ y 1 -GiuMW-GW) 



(24) 



for all a € [0, °°) and for some y > 1 . Moreover, the right end of the interval above is 
given by 



Um^(a) = |b"-G(i4) 



(25) 



Since Km(oc) is continuously increasing, it follows from (24i that there exists a* G 
[Q,°°) such that 



-lb" -G(ul)\\ 2 < p 2 \\yii-G(uZ)\\ 2 = da*) < -G(«£) 



(26) 



Note that any a„, such that a* < a m will therefore satisfy fc„(a*) < fCm(a m ) which, 
from ( po) implies |9} as required. Computationally, we can determine such a m by 
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constructing aj„ — s- °° as j — »• oo, Let us consider, for example, a,i +1 — 2^ +x aL where 
> is an initial guess for a m . We claim that there exists J < °° such that 



If no such J exists then 



p 2 |b"-G(4)H 2 <^«) (27) 
^{aL)<p 2 \\y^-G{ul)\\ 2 (28) 



for all j £ N. In particular, for sufficient large j, from ( 25 i we find 



lb T '-G(4)ll 2 <p 2 lb , '-G(«3[)ll 2 (29) 

which contradicts the hypothesis of p < 1. We define a m = a ] m and the update of the 
regularizing LM scheme 

«*H = ul+CDG*(ul)[DG{ul)CDG*(ul) + a m r]- x [y*-G{ul)] (30) 

Note that the computation of a m requires the evaluation of K„ (a) which from ( 23 i in- 
volves the inversion of [DG(u%,) C DG* (ui) + a m r] - 1 . However, DG(u} n ) C DG* (u? n ) 
has to be assembled only once per iteration of the scheme (see the update equation 



(30 1). The cost of inverting [DG(uH) C DG*(uJj,) + a m F] _1 for different a m 's is neg- 
ligible for the application under consideration due the small dimensionality of the 
observation space. Therefore, the cost of computing a m that satisfies (|9| is negligible 
compared to the cost of evaluating G(um) and assembling DG(4) C DG* (um) which 
both, in turn, constitute the main computational cost per iteration of the proposed 
implementation of the regularizing LM scheme. With the aforementioned consider- 
ations we propose a computationally efficient implementation of the regularization 
LM scheme. 



Algorithm 2 (Regularizing LM Scheme (implementable version)) Let uq GX be 

an initial guess. Choose parameters p £ (0, 1) and T > 1 /p. For each n = 1 , . . ., 

(1) Solution to the forward model. Given u\ evaluate the forward operator G(um). 

(2) Stopping rule (Discrepancy Principle). If 

\\r l l 2 {y^-G{ul))\\y<x7] (31) 

stop. Output: u? n . 

(3) Compute the sensitivity matrices DG(um), its adjoint operator DGium)* and as- 
semble matrix DG(um) C DG* («,?,). Let > and cc/^ 1 — 2 ; a,i +1 . Let J be such 
that 

p 2 ||y"-G(^)|| 2 <^«) 
= a 2 \\r l l 2 \DG{ul) CDG* {ul) + a J m r] 1 \f - G(ul)] \\ 2 (32) 

Update. Define 

ul +l =ul+CDG*{ul)[DG{ul)CDG*(ul) + a J m r]-\y*~G{ul)] (33) 
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3.3 The standard approach for history matching 

As we indicated in SectionfT] one of the most standard approaches for history match- 
ing consist of minimizing l2l). For analogy with our proposed implementation for 
solving ([TJ, we based the following discussion on the the application of the Levenberg- 
Marquardt algorithm used in [16,21 22 , 19] for the minimization of |2]) in the stan- 
dard approach. The aforementioned method consist of computing the sequence u m+ \ = 
u m +Au where the step Au satisfies 



DG* (u m )r- l DG(u m ) +C" 1 + KfT x Au=DG*{u n )r- l [y T] -G{u m ))-C-\u 



-lr..J7 nt.. \\ ^-l. 



for some A m > 0. The proposed update (34i is a scaled version of the standard LM 
algorithm for the solution of well -posed optimization problems [18. Chapter 10]. For 
the history matching applications of [16,21 22, 19|, the suggested selection of A m is 
the following. The initial Xq is chosen between -J 'J(uq) /Nj and J(uq)/N 1 i where Nj 
is the dimension of the observation space. For m > 0, Am+i is chosen according to 

1 _ / Am/10 if J{u m+ \) < J{u m ) 

Am+1 -\ 10A,,, if J{u m+l )>J{u m ) W 

In addition, the stopping criteria for the LM technique of [ 16,21 .22 . 19 1 is based on 
the following two stopping criteria 

\J(um+l) ~J{u 



J(u m -\-\) 
\\ u m+l — u m\\x 



< £o, (36) 
< ei (37) 



\\ u m+i\\: 

In order to understand the standard LM approach for minimizing (|2]), note that ( 34 1 
can be derived from the Euler-Lagrange equations for the minimization of 

Q m (Au) = l -\\r- l ' 2 (y n -G{u m ) -DG{u m )Au)\\ 2 

+ 1 -\\C-^ 2 (Au~{u m -u))\\l+ l -^ n \\C- l/2 Au\\ 2 x (38) 

In other words, Au m = argmin ve x Q' n (v), and so each iteration step of the LM method 
of 02 ,21 ,22,19] is the solution of a least-squares Tikhonov-type problem on the 
linearized inverse problem. Note that the choice A,,, = in (38 1 suppresses the extra 



regularization term in the right hand side of ( |38| . Indeed, the initial motivation of the 
LM scheme used in [ 16,21 ,22, 19] was to alleviate the lack of stability of the Gauss- 
Newton (GN) method of IT6l and lfl9l section 8.4.2] which corresponds to A m = in 
(|34|> ((38). The LM scheme of lfT6ll2T1l22l[T9l for the minimization of ^ is an efficient 
strategy provided that the minimization of / is a well-posed problem. However, as 
we indicated before, the regularization term R(u) in (|2| may be insufficient for some 
choices of C. For some choices of C, in the following section we present numerical 
experiments demonstrating that the selection of A in ( |3~5| ) and the stopping criteria 
of <j36j-(j37j may lead to both lack of stability and fidelity in the computation of 
estimates of geologic parameters. 



12 



Marco A. Iglesias, Clint Dawson. 



3.4 Computational similarities between the standard and the proposed approach 



We emphasize that the regularizing LM scheme presented in Section [3] is designed 
to compute stable approximation to a minimizer of <P defined in ([T]). In contrast, the 
standard approach discussed in the preceding section is based on minimizing Q by 
means of a standard optimization algorithm. Therefore, the two approaches aim at 
solving two substantially different problems. Nonetheless, there are some computa- 
tional similarities between the aforementioned approaches as we now discuss. Let us 
recall that each iteration step for minimizing (|2]i in the standard approach is given by 
d34| which, from the lemmas in |fT9"l Section 7.4], is equivalent to 



Au = CDG*(u m ) DG(u m )CDG*(u m ) + {l+X m )r y n -G(u m ) 



-1 r 



+ 



1+A, 



-DG(u m )(u m -u) 



+ 



1+A, 



-{u m -u) 



(39) 



On the other hand, the mth step computed with the regularizing LM scheme for ap- 
proximating the minimizer of ([TJ is given by 

Aul = CDG* {uZ)[DG(ul)CDG* (n£) + a m r] ~ l \f - G{ul)] (40) 



The substantial similarities between expression (39 1 and (40 1 are evident although 



they converge to different functions. Notice that, at the discretization level, the com- 
putation of the matrices CDG*(u) and DG(u)C DG* (u) as well as the evaluation of 
G(u) are needed for both approaches (obviously evaluated at different m's). In addi- 
tion, note that the terms 1/(1 + X m )DG(u m )(u m — u) and 1/(1 + X m )DG(u m )(u„ 

u) are not required in ( 40 1. It is therefore clear that the main routines and codes 



used for computing (39 in the standard approach can be used for implementing 



the regularizing LM scheme step (40 1. In fact, a routine that assembles CDG*(u) 



and DG(u)CDG*(u) as well as the routine that evaluates G(u) are sufficient for a 
straightforward implementation of the regularizing LM scheme Algorithm [2] The 
aforementioned computational similarities open the possibility to study the history 
matching problem in the sense presented in this paper by using available implemen- 
tations for the standard approach. Moreover, provided that the same implementation 
for CDG* (m), DG(u) CDG* (u) and G(u) are used for both approaches, from the pre- 
vious discussion it follows that the two approaches have the same computational cost 
per iteration. For the results presented in the subsequent section, the operator DG and 
DG* are computed as described in Section 9.7 of |fT9l . 



4 Numerical Results 



In this section we present numerical experiments to show the capabilities of the regu- 
larizing LM scheme for estimating the log-permeability in the oil-water model of the 
Appendix. 
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4. 1 Experimental setting 

We consider a synthetic experiment where the reservoir domain is £2 = [0,L] x [Q,L] 
and the prior knowledge of the subsurface is given in terms of a prior u = 500 md 
(constant in £2) and a covariance operator 

C = K- l C (41) 

where Co is a spherical covariance function Q 

( i 3 \\M g ( Zi -z 2 )\\ , l Hj^fa-ajHj i f iiM r 7 , 7^11^/7 

[ if ||M e (zi-z 2 )|| >a 

with z, = (x;,y,). In the previous expression, Mg is a rotation matrix along the direc- 



tion of maximum continuity with range denoted by a. Covariance functions like (42 1 



are common in modeling geologic properties of reservoirs 13] . The tunable parame- 



ter K in (42 1 will enable us to study the performance of the proposed approach with 
respect to different choices of the prior covariance parameterized in terms of K. 

We consider K = 1 in ( |4"Tj i to be the "correct" covariance in the sense that the 
true (or reference) log-permeability is a Gaussian field with mean u and covariance 
C = Co- In other words, K = 1 corresponds to the best case scenario where with our 
prior knowledge is consistent with the truth. In Figure [T] (left) we display the true 
permeability u\ sampled from the aforementioned distribution. We now consider a 
water flood described with the model presented in the Appendix. Nine production 
wells P\,...,Pg and four injection wells I\ , . . . J4 are considered in the configuration 
displayed in Figure [T] (right). Relevant data of the reservoir model is displayed in 
Table [T] We use the true log-permeability field of Fi gure [T] (l eft) to generate syn- 



thetic data as we now describe. First, the PDE system (46 i-|47]i is solved for u 
the resulting pressures and saturations are used in the expression for the measure- 
ment functional (53 1-(56 1 to find G(uf). Finally, synthetic data is generated by adding 



Gaussian random noise \ ~ N(Q,r). More precisely, we define y 1 ! = G{u*) + £. We 
consider a diagonal error measurement covariance F with diagonal elements denoted 
by of. The values of d, associated to measurements of bottom hole pressure consist 
of some percentage (defined below) of the nominal value of the corresponding mea- 
sured variable. In order to avoid zero values for the (7, 's associated to measurements 
of water rates, for either water and oil rate measurements, the corresponding a,- is a 
percentage of the nominal value of the total flow rate (which is the well constraint). 
The aforementioned percentage is the same for both measurements of pressure and 
flow rates. The noise level is defined by 

^llr-^V-G^))^ (43) 



4.2 Performance of the LM scheme with respect to the observational noise level 



In this subsection we investigate the accuracy of the estimate obtained with the 
regularizing LM scheme as a function of the noise level rj with 77 — > 0. According to 
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Table 1 Reservoir model description 



Variable Value | Variable Value 



L[m 3 ] 


2x 10 3 


a w 




0.3 


c [Pa" 1 ] 


0.0 


a a 




0.9 


v Q [Pa s] 


10- 2 


b pi 
r bh 


[Pa] 


2.7 x 10 7 


T [years] 


5 


h l 


[m 3 /day] 


2.6 x 10 3 


a Po [Pa] 


2.5 x 10 7 


Siw 




0.2 




0.2 


$ro 




0.2 


V w [Pa s] 


5 x 10~ 4 









a Constant in 13. b Constant in [0,7]- 



Truth well locations 




x (km) 



Fig. 1 Left: True log-permeability [logm 2 ]. Right: Well configuration. 



Theorem[T] ifl converges to a minimizer of ([TJ as rj — » (see Remark|2|. Although 
this converged solution may not necessarily be the truth u* (due to possible non- 
uniqueness), the solutions may arguably reflect the main spatial features of the truth. 
We therefore consider the accuracy of the estimates in terms of their relative error 
with respect to the true log-permeability w . 

Five sets of synthetic data associated to different noise levels {Tjy}y =1 are gener- 
ated with the procedure previously described. For each set, a different percentage of 
the nominal value of the measured values is selected. The resulting sets of synthetic 
data {y^}^! provide noise levels (defined by that correspond to some fractions 
of the norm of the corresponding measurements | |F _1 / 2 -y1j j G {1,...,5}. More 
precisely, we have 

rjj^fjWr- 1 ^]] (44) 

with /i = 5 x 1(T 2 , f 2 = 1(T 2 , / 3 = 5 x 1(T 3 , f 4 = and / 5 =5x 1(T 4 

For this set of experiments, the parameters for the regularizing LM scheme are se- 
lected as T = 1.2 and p = 0.83. Further choices of p and T are investigated in subsec- 
tion 4.3 In addition, we consider K = 1 in (41 i.The performance of the regularizing 
LM scheme for each of the five sets of synthetic data corresponding to different noise 
levels is presented in Figure^ The data misfit | |F _1/ ' 2 (y T ' — G(u^))\ |y is displayed in 
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2 1 1 1 1 1 1 1 1 0.65 ' ' ' ' ' ' ' 1 

□ 5 10 15 20 25 30 35 40 5 10 15 20 25 30 35 40 

iteration number iteration number 

Fig. 2 Performance of the regularizing LM scheme with respect to the noise level . Right: data misfit. Left: 
relative error with respect to the truth 



Figure [5] (left) and the relative error with respect to the truth ||m„,+i — u m \ |x/|| M ;«+i ||x 
is shown in Figure [2] (right). The stability in the computation of the numerical solu- 
tions is reflected in the decrease of the relative error with respect to the truth. Note 
that, as the noise level decreases, the accuracy with respect to the relative error in- 
creases. The dependence of the accuracy on the noise level can be visually appreci- 
ated from the log-permeability estimates presented in Figure [3] For smaller noise in 
the observations, the regularizing LM scheme seems to provide stable and accurate 
estimates of the geologic properties. 

By construction, the weight F~'/ 2 (in the data misfit) depends inversely on the 
error y* 1 —G{u*). Therefore, even though the five experiments have the same initial 
guess mo = u, the initial value of the data misfit is larger for smaller noise levels. Fur- 
thermore, since the error that defines the noise level d4) is also weig htedbyr-^fhe 
actual value in (04) is approximately similar for all the experiments. The difference, 



however, is in the corresponding fraction fj of the norm in ( 44 1 which is used in the 
label of Figure [2] 



4.3 Parameters T and p 

The parameters p e (0,1) and T > 1/p are the tunable parameters in the regularizing 
LM scheme. In this section we present the numerical performance of the LM scheme 
for different choices of these parameters. We recall from Section [3] that if p « 1 
(p < 1), we then may choose T ~ 1 (t > 1/p). Then the regularizing LM scheme 
terminates when the estimate u} n produces a data misfit H-T^ 1 / 2 ^ 1 ' —G(u}^\\y ~ rj. 
Small values of p imply larger values of T which may lead to estimates that provide 
a poor fit to the production data. It is therefore important to study the potential lack 
of accuracy due to the choices of p and T. 

We consider the same experimental setting as before for only one fixed set of 
synthetic data with 1 % of observational noise level. We consider several choices of 
p, with the corresponding T defined by x — l/(p — 10 ^). The performance of the 
LM scheme for these choices of parameters is presented in Figure HJ and Figure [5] 
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Fig. 3 Log-permeability estimates obtained with the regularizing LM scheme for different noise levels 
[(logm 2 )] 




Although reasonable estimates were obtained for all these choices of p, from Figure 
|4]we observe that more accurate estimates, in terms of the relative error with respect 
to the truth, are obtained when p is indeed close to one. However, it is important 
to remark that an increase in the computational cost is associated with the improved 
accuracy for p ss 1. These numerical experiments suggest that optimal choices in 
terms of computational efficiency and accuracy are obtained for p £ [0.8,0.9]. 
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1 
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x(km) 

p=0.8 



f 



p=0.95 



x(km) 

p=0.75 



ill 



p=0.9 



■ 



x(km) 

p=0.65 



Fig. 5 Log-permeability estimates obtained with the regularizing LM scheme for different parameter p. 
[(logm 2 )] 



4.4 Performance of the LM scheme with respect to the prior covariance. 

Recall that in the previous experiments we have chosen K = 1 in ( pT) which corre- 
sponds to the best-case scenario where the true log-permeability is consistent with the 
prior knowledge. In this subsection we investigate the performance of the regulariz- 
ing LM scheme with respect to different choices of the prior covariance. In particular, 



we consider the case where the prior covariance is parameterized by (41 1, and we are 
interested in the performance of the LM scheme with respect to K. The values of the 
parameters z and p in the regularizing LM scheme are as described in the subsection 



4.2 and the synthetic data as the same as in subsection 4.3 



In Figure [6] (left) we present the data misfit associated to the regularizing LM 



scheme for some choices of K with K < 1 in the prior covariance (41 1. The hor- 
izontal line indicates the value of T7] used in the stoping criterion. We recall that 
the regularizing LM scheme is stopped after the estimate produces a data misfit be- 
low the aforementioned value. In Figure [6] (right) we display the relative error of the 
estimates with respect to the truth. In Figure [7] we present the performance of the 
regularizing LM scheme for K > 10 in (41) . The log-permeability estimates for all 
k's are displayed in Figure [8] It is clear that the regularizing LM scheme produce 
similar estimates regardless the value of K in the covariance expression ( |4T) . For the 
present experiments, the estimates reach the stopping criterion after approximately 



15 iterations. In Figure [9] and Figure 10 we display the model predictions obtained 



by simulating the water flood with the estimates of log-permeability produced with 
the regularizing LM scheme. As we expect from the similarities in all the estimates 
(Figure |8) obtained for all K's considered here, the associated model predictions are 
all almost identical. 
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2 4 6 B 10 12 14 16 5 10 15 



iteration number iteration number 

Fig. 6 Regularizing LM scheme for history matching. Right: data misfit. Left: relative error with respect 
to the truth 




2 4 6 8 10 12 14 16 ' 2 4 6 8 10 12 14 

iteration number iteration number 

Fig. 7 Regularizing LM scheme for history matching. Right: data misfit. Left: relative error with respect 
to the truth 



4.5 Comparison with the standard approach 



We consider the same set of synthetic data yt , measurement error covariance F, prior 



mean u and C (for the same jc's) used in the experiment of subsection 4.4 In this case, 
however, we find estimates of the log-permeability by means of the standard approach 
of [16 .21 22 19 1 described in subsection 3.3 Note that with the choice of C given 



by (41 1-(|42")>, the objective functional that is minimized in the standard approach (j2j) 



becomes 



/(«) = -||r- 1 /2(y«l_ G(M)) ||2 +)c . | (C5 -1/2 (B 



(45) 



Therefore, K in (41 1 controls the relative size the the prior term with respect to the 



data misfit. In Figure 1 1 we report the performance of the experiments for K < 1 . 



The right panel of Figure 1 1 shows the relative error with respect to the truth of 
the estimate log -permeability field. In Figure [TT] (left) we present the associated log- 
objective functional ( |4"5j ). As the number of iteration increases, the method produces 
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Fig. 8 Log-permeability estimates obtained with the regularizing LM scheme for different fc's in 
(refeq:4.19) [(logm 2 )] 



estimates that decreases the objective functional /. However, due to the lack of sta- 
bility in the computation of ( 39 1, the error with respect to the truth increases after 
a certain number of iterations. Note that, even when the estimate is computed with 
(k = 1) the same covariance used for the generation of the truth, the corresponding 
error starts increasing after 5 iterations of the method. Additionally, Figure 1 1 reveals 
the potential failure of the stopping criteria (|36|-([37]i in the standard approach of |[T6l 
l2Tll22l[T9l . More precisely, due to the ill-posedness of the inverse problem, a decrease 
of the objective functional (j2j) may not be associated with a controlled change in the 
corresponding estimate u. Therefore, the choice of A based on (35 1 may still lead to 
large values of the estimate for which ( 37 1 may not be satisfied. 

We now consider the minimization of d45j for K > 10. The effect of larger reg- 
ularization in (45 i is observed in Figure [l2](right) where, after some number of it- 
erations, the error is indeed stabilized. In Figure 12 (left) we show the associated 
objective functional. In this case, the stopping criteria |36]l and ( [37] i are both satis- 
fied. However, for larger fc's less accurate estimates are obtained. The estimates of 
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time (years) time (years) time (years) 



Fig. 9 Water rates [bbl/day]. From left to right: Wells P%, P4 and P5. Top: Experiments for K < 1. Bottom: 
Experiments for K > 10 




time (years) time (years) time (years) 



Fig. 10 Bottom hole pressure [Pa]. From left to right: Wells I2, h and I4. Top: Experiments for K < 1. 
Bottom: Experiments for K > 10 



the log-permeability obtained for all it's after 35 iterations of the standard method 
are displayed in Figure [13] For small K, the lack of stability in the computations is 
reflected in very large values of the log-permeability fields which is consistent with 
the results reported in 11611191 . On the other hand, for large K, the lack of fidelity of 
the corresponding estimates can be visually appreciated for K > 10 2 . From Figure 



1 1 (right) and Figure 12 (right) we conclude that the correct choice of C (i.e. K = 1) 



does not lead to the optimal estimate in terms of the error with respect to the truth. 
In fact, from all the experiments, K — 10 2 provides the minimal error with respect to 
the truth. Similar to the previous set of experiments, in Figure [14] and Figure [TB] we 
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iteration number iteration number 

Fig. 11 Performance of the standard approach for history matching (if < 1). Left: Objective functional 
|2j. Right: Relative error (left hand side {57}) 



show model predictions during the 10 years of history matching and the prediction 
time of 5 years. Note that for small K (k < 1), all the estimates provide a good data 
match even though the quality of the corresponding geologic properties (see Figure 



13 i is severely degraded for small K. In contrast, the lack of fidelity for larger values 



of K corresponds to poor estimates of the data match as we expected. 



In contrast to the standard approach (see Figure 13 I, even for small K in (41 1 the 



regularizing LM scheme produce stable estimates of the the true log-permeability 
(see Figure [8]). In addition, the geological constraint of C is enforced with the regu- 
larizing LM scheme. However, by producing stable estimates of the minimization of 
([l} the LM scheme avoids the potential lack of fidelity of the standard approach due 
to the potential overestimation of the prior term. Moreover, as we indicated earlier, 
the computational cost per iteration of our implementation of the regularizing LM 
scheme is equivalent to the cost per iteration of the standard method of 116 2T1I221 
[T9l for minimizing (Kb. From Figure [l2j and we observe that for K > 10 2 , the con- 
vergence for both approaches is achieved after 15 iterations. However, the accuracy in 
terms of the relative error seems to be outperformed by the regularizing LM scheme 
for larger values of K. On the other hand, for the case with K < 1 , the convergence 
of the standard approach is not achieved due to the lack of stability reflected in the 
increase in the relative error. 



5 Conclusions 

While the main contribution of this paper is the implementation of the regularizing 
LM scheme for history matching, our general aim is to promote further investigations 
of implementation based on well established theories for approximating stable solu- 
tions of history matching problems posed as the minimization of ([TJ. Although the 
discussions and experiments of this paper are based on the LM method, the funda- 
mental ideas can be applied to other techniques. In particular, in the iterative regular- 
ization literature, there are analogous gradient based (e.g. Landweber iteration, steep- 
est descent) and quasi-Newton methods (BFGS, conjugate gradients) whose aim is 
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Fig. 12 Performance of the standard approach for history matching (k > 10). Left: Objective functional 
{2). Right: Relative error (left hand side )37) ) 



Truth 



k=5x 10 



K=1(T 




Fig. 13 Log-permeability estimates obtained with different k's in |^J (i.e. the standard approach for history 
matching) [(logm 2 )] 
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time (years) time (years) time (years) 



Fig. 14 Water rates [bbl/day]. From left to right: Wells Pi, P4 and P5. Top: Experiments for K < 1. Bottom: 
Experiments for K > 10 




time (years) time (years) time (years) 



Fig. 15 Bottom hole pressure [Pa]. From left to right: Wells I2, h and I4. Top: Experiments for K < 1. 
Bottom: Experiments for K > 10 

to solve nonlinear inverse ill-posed problems such as the one presented here. More- 
over, similarly to the LM scheme, those gradient-based and quasi-Newton methods 
also share similarities with the corresponding ones for solving (|2]i in the standard ap- 
proach. Those similarities may lead to straightforward implementations of iterative 
regularization techniques when the standard technologies are already available. 

Conducting history matching by minimizing Q has been often motivated from 
the Bayesian formulation for data assimilation. Under Gaussian assumptions, the 
minimizer of also called maximum a posteriori (MAP) estimate, maximizes the 
conditional posterior probability measure of the unknown given the observed data 
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|[T9l . While our proposed history matching approach based on iterative regulariza- 
tion techniques is entirely deterministic, there is a potential use of these techniques 
within the context of Bayesian data assimilation for uncertainty quantification. This 
conjecture follows from the fact that some standard techniques that approximate the 
posterior distribution of the Bayesian framework are constructed by randomizing the 
solution to deterministic problems [ 10 1 . On the other hand, iterative regularization 
techniques can also be used in the context of facies identification as suggested in 
ifTTI where a geometric -based iterative regularization approach was applied for the 
estimation of geologic facies given data from an oil-water reservoir model similar 
to the one considered here. Iterative regularization provides then a broad spectrum 
of techniques that can be potentially used in the estimation of geologic properties in 
reservoir models. 



6 Appendix: The forward operator 

We recall that G is the forward operator that maps the geologic parameters to the 
production data. We briefly describe the forward operator that we use for the nu- 
merical experiments presented in Section 3.3 and Section[4] We consider simplified 



two-dimensional models typical for testing history matching algorithms. The domain 
of the reservoir is denoted by D; the absolute permeability and porosity are denoted 
by K and respectively. The interval [0, T] (T > 0) is the time interval of interest for 
the flow simulation. For simplicity we assume that the only unknown parameter is 
u = \ogK. Nevertheless, all the techniques and implementations that we describe in 
subsequent sections can be extended to include additional parameters (e.g. porosity). 

We consider an incompressible oil-water reservoir model initially saturated with 
oil and irreducible water, the water and oil phase are indexed by j3 = w and j3 = o, 
respectively. We are interested in a waterflood process where water is injected at Ni 
injection wells located at {jc[}i=i ■ Water and oil are produced at Np production wells 
located at{jCp}^j. Additionally, we assume that injection wells are operated under 
prescribed rates {q\{t) }._£j while production wells are constrained to the total flow 

rate {qp(t)},Zi- The pressure p(x,t) and the saturation s(x,t) ((x,t) E D x [0, T]) are 
the state variables. From standard arguments it can be shown that (s,p) is the solution 
to the following system U 

N, N w 

- V • A (s)e u Vp = £ q\8{x - x[) + £ q' P 8(x - x' P ) (46) 
1=1 1=1 

- - V • X w (syVp = t q^Six-Xj) + £ ^q P 8(x-x l P ) (47) 

in D x (0, T], where 8(x — Xp) and 8(x — x\) are the (possibly mollified) Dirac deltas. 
In (|46]i-(|47)>, ^(s) and A (s) denote the water and total mobility defined by 

K{s) = KM = lUsl + 
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where k r y(s) and denote the relative permeability and the viscosity of the y-phase 
fluid, respectively. Furthermore, we assume that 



s - su 



1-5 



iw °or 



1 — S — S n 



1-s 



iw ^or 



(49) 



where a W) a € (0, 1], Si w is the irreducible water saturation and s or is the residual oil 
saturation. We additionally prescribe initial conditions for pressure and water satura- 
tion 

p = P0, s = s inDx{0} (50) 
For simplicity, no-flow boundary conditions are prescribed on the reservoir boundary 



-e"X{s)Vp-n = ondDx(0,T] 
~e l % v (s)Vp-n = ondDx(Q,T] 



(51) 
(52) 



Let us assume that there are Nm measurement times denoted as before We 
assume measurements of bottom-hole pressure are collected at the injection wells at 
{ { >'}n=i - This, according to Peacemen well-model [1 1 is defined by 



*•'<'"> = «S> +P, -< U) 



(53) 



for I — 1 , . . . , Nj and n = 1 , . . . ,Nm- Analogously, we consider measurements of water 
and oil rates at the production wells 



X(s(x' p ,t n )) ' " ' X(s(xf P ,t„)) 



q l P (tn) (54) 



for I = 1, . . . , Np and n = 1, . . . ,N&i- In ( |54| , A„ = A - A vr . Let us define the 2Np +N[- 
dimensional vector 

M n (p,s) = (My(p,s),. . . ,M^{p,s),M l n ^{p,s),. . . ,M^{p,slMl A M, ■ ■ ■ ,M»"> p °(p,s)) 

(55) 

that contains the number of measurements from wells at a given time. The total num- 
ber of measurements is N — [2Np +Nj]Nm and the forward map G:X-> M. N is then 
given by expression 



G{u) = (M\ (p,s),... ,M Nm (p,s)) 



(56) 



which comprises the production data obtained from production and injection wells at 
the measurement times. 
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